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ABSTRACT 

We invert directly the redshift - luminosity distribution of observed long Swift GRBs 
to obtain their rate and luminosity function. Our best fit rate is described by a broken 
power law that rises like (l+z)^'^-"6 forO < z < 3 and decrease like (l+z)^^''-io 
for z > 3. The local rate is po — ^■•^-o'rlGpc^'^yr^^]. The luminosity function is 
well described by a broken power law with a break at L* ~ j^q''^- ''^"■^ [erg /sec] and 
with indices a = 0.2+°;^ and /3 = lA^oi- The recently detected GRB 090423, with 
redshift w 8, fits nicely into the model's prediction, verifying that we are allowed to 
extend our results to high redshifts. While there is a possible agreement with the star 
formation rate (SFR) for z < 3, the high redshift slope is shallower than the steep 
decline in the SFR for 4 < z. However we cannot rule out a GRB rate that follows 
one of the recent SFR models. 

1 INTRODUCTION 

Gamma-ray bursts (GRBs) are short and intense pulses of soft 7-rays. In this work we study their 
luminosity function and their cosmic rate. These functions are essential to understand the nature of 
GRBs and to determine their progenitors. They may also shed light on the still mysterious physics 
of the central engine. 

The luminosity function and the cosmic GRB rate are observationally entangled as the ob- 
served rate is a convolution of the luminosity function with the cosmic rate. For this reason, almost 
every work before this paper have made an a-priori assumptions on at least one of these functions. 
At first - having no motivation to believe otherwise - the rate was assumed to be constant. The 
simplest form for the luminosity function was a standard candle i.e. a constant luminosity. These 
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early studies used the measured (V^/Knax) value to find the typical luminosity (Mao and Paczyn- 



ski|1992| |Piran|1992| [Fenimore et al.|1993^ |Ulmer and Wijers|1995[ |Ulmer et al.fT995] ) and corre- 



spondingly a maximal distance from which GRBs were observed. Later, using the flux distribution 



(logN — logP relation), [Cohen and Piran| ( fT995| ) and |Loredo and Wasserman| ( |1998| ) showed how 
relaxing the standard candle assumption allows a corresponding relaxation of the constant rate 
assumption. 

[Paczyhski ( |1998 1 noticed that the hosts of the bursts are in star forming regions and suggested 
that GRBs follow the SFR (see also Totani 1997 [Wijers et al.|1998al ). The detection of a supernova 
associated with GRB980425 ( |Galama et"aL]|1998| ), strengthened the expectation that the GRB 
rate should follow the SFR. Using this proportionality, numerous studies examined the typical 



luminosity assuming at first standard candles (e.g. Wijers et al. 1998a, Totani 1999) and later 
more elaborate shapes of the luminosity function (e.g. |Schmidt|1999[|2001b|a[ |Guetta et al.|2005 
Firmani et al.|2004t[Guetta and Piran|2005] ). 



Swift that discovers routinely GRBs afterglows detected GRBs from higher redshifts than was 



previously possible sparked renewed interest in the GRB redshift distribution (e.g. Berger et al. 



2005 , Natarajan et al.||2005 , Bromm and Loeb|2006 


Jakobsson et al.||2006 , Le and Dermer|2007 


Yiiksel and Kistler |2007[ Salvaterra and Chincarini 


2007, Liang et al. |2007t Chary et al.|2007|). 



More and more signs appeared suggesting that the rate of GRBs does not simply follow the global 
star formation rate, as was believed earlier. [Firmani et aL] ( |2005[ ); |Le Floc'h et al.| ( |2006[ ); [Daigne 



et al. (2006); Le and Dermer (2007); Guetta and Piran (2007) conclude that the GRB rate differs 



from the SFR. In particular we observe more high redshift bursts than what is expected for a rate 



that follows the SFR. Yiiksel et al. (2008) uses the high luminosity subsample of bursts to obtain 
the GRB rate without assuming any luminosity function. They find that the GRB rate at high 
redshifts (z > 4) is higher on than Hopkins and Beacom| (2006) SFR. Alternatively luminosity 
evolution, has been suggested by some authors: (e.g. [Lloyd-Ronning et aL][2002^ [Firmani et al. 



2004|[Matsubayashi et al.|2006t[Kocevski and Liang|2006[ ). [Salvaterra et~aL]p008| ) used a sample 



of long bursts with redshift and found evidence for luminosity evolution. Other papers, including 
this one, find self consistency without a luminosity evolution. 

We introduce here a new method for determining the rate and the luminosity, by inverting the 
observations without making any assumptions on the functional form of the luminosity and rate 
functions. This direct inversion allows us to use most of the available redshift data and to obtain 
robust estimation of both functions. 

In ^we describe the sub-sample of bursts we analyze and its advantages over other samples. 
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In ^ we introduce a new formalism allowing us to directly invert the observed distribution and 
obtain the intrinsic luminosity function and rate. In ^|4] we present the results obtained using this 
method and the confidence ranges associated with them. In ^we address the question whether 
the GRB rate is consistent with the SFR. We discuss a few implications of the results in ^ 



2 THE SAMPLE 

We consider long bursts (tgo ^ 2 sec) detected by Swift from the beginning of its operation until 
burst 090726 with a measured peak-flux and a measured redshift. Generally GRBs redshifts are 
obtained from the optical afterglow spectrum using absorption lines or photometry, or from the 
spectrum of the host galaxy using emission lines. However in the global sample we find differ- 
ent redshift distributions for the different detection methods, namely: Absorption, Emission and 
Photometry (see Figure [T]). For redshifts determined using the hosts' emission lines, we do not 
detect high-redshift events, whereas absorption lines redshifts extend over the entire range of red- 
shifts. Furthermore, emission lines are more susceptible to a selection effect known as the 'redshift 
desert' in the range 1.1 < z < 2.l[^ ( [Fiore et al.|2007 , see also Coward|2008 1. This is in line with 



the fact that we have also found that the probability to measure the redshift using emission lines 
strongly depends on the gamma ray flux, favoring high fluxes This effect is mild for absorption 
lines and photometry (see appendix B). 

To obtain a more uniform sample we consider therefore only bursts whose redshift was mea- 
sured using the afterglow. For each burst we calculate the isotropic equivalent peak-luminosity: 
Liso (see Appendix A. for details) using the peak-flux and redshift. We use standard AC DM cos- 
mology with h = 0.7, fi^ = 0.27, Qa = 0.73. 



3 A DIRECT ESTIMATE OF THE LUMINOSITY FUNCTION AND THE RATE 
3.1 Assumptions 

We assume that the luminosity function is redshift independent. In this case the number of bursts 
at a given redshift and with a given luminosity is the product of the luminosity function, 0(L), that 
depends only on the luminosity L and the GRB rate, Rgrb{z), that depends only on the redshift 
z. This common assumption is reasonable since a-priori there is no competing reason why the 



1 

2 

3 



All data were taken from the Swift information page http : / / swift.gsfc.nasa.gov/docs/ swift/archive/ grbJ^able/ 
Although a redshift determination through absorption lines is also difficult in the range 1.5 < z < 2.1 
A possible explanation is that a brighter burst can be more easily localized, making follow-up possible. 
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Figure 1. The redshift distribution for the different methods: Absorption lines, Emission lines and Photometry 



luminosity function should depend on the redshift. We test later the validity of this assumption and 
show that it is accepted with a high statistical significance. 

The isotropic peak luminosity function, is defined traditionally as the fraction of GRBs 
with isotropic equivalent luminosities in the interval log L and log L + d log L. The rate, Rgrb{z), 
is defined as the co-moving space density of GRBs in the interval z io z + dz. The distribution 
density: n(L, z) is given by: 

n{L, z)dlogLdz = (f){L) ■ R{z)dlogLdz , (1) 

where 

Rgrb{z) dV{z) 



R{z) 



(2) 



(1 + z) dz ' 

is the differential co-moving rate of bursts at a redshift z, dV{z)/dz is the derivative of the volume 
element and the factor (1 + z)^^ reflects the cosmological time dilation. 



3.2 Formalism 

We consider now a direct method to invert the observed L — z distribution and obtain the functions 
and Rgrb- To do so we approximate (f){L) and R{z) as step functions whose range is divided to 
bins with a constant value within each bin. These functions can be expressed as a sum of Heaviside 
functions. 
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First, we divide both luminosity and redshift intervals into bins. We denote 

0i = ^ L < Li+i) , (3) 



[zj+i - Zj) J,^ (1 + z) dz 

We also define the weights factors 



Wij = / 6{L, z)dlogLdz , (5) 

J Li J Zj 

as the probability for detecting a burst with a measured redshift z and luminosity L where 9{L, z) = 
9z{p{L, z)) is the probability to detect and measure redshift for a burst with a luminosity L at a 
redshift z (see appendix B). We denote the observed number of events per bin as: 

= N{U ^L< , (6) 

= N{z.j <:z< Zj+i) , (7) 

Nij = N{Li^L< Zj ^z< Zj+i) , (8) 

iV^J^iV,,. (9) 

Next, we determine (pi and Rj and the error estimates, using the maximum-likelihood formal- 
ism. We define M the (logarithmic) likelihood of the model given the observations as: 

M = Y1 H(f>iRjWij] - N ln[^ (piRjWij] . (10) 

At the maximum all partial derivatives of M with respect to (pi and with respect to Rj vanish, 
leading to 

Ni Y.i'j'(pi'RfWi'j' 

(l>i = , (11) 

and 

Nj Y.i>ji (pi'Rj'Wi'ji 
J2i (t^iWij N 

Notice that the second term on the right hand side of each of the equations 1 1 and 12 is the same 
normalization factor. We have a set of non-linear equations with as many equations as variables. 
We solve numerically these non linear equations using successive iterations until convergence. A- 
priori it is not clear whether (p and R are uniquely determined and whether there is a solution at all. 
However, we find good convergence. We have examined a large set (10^) of initial guesses where 
each component was randomly drawn from a uniform distribution. We found a rapid convergence 
to a unique solution for all initial guesses, all reaching the requested accuracy of 10^^ with less 
than 25 iterations. We thus conclude that the existence of other stable solutions is very unlikely. 
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3.3 Error Estimates 

We approximate the error as the value for which M deviate by -1 from it's maximum: (i.e. the 
likelihood is smaller by a factor e). This reflects an Icr error for a normal-distribution. 

-1 = N, In 1 + - iVln 1 + , (13) 

(pi N (pi 

- 1 = AT, ln(l + ^) - NHl + • (14) 

The two solutions, i.e. the positive one and the negative one, give an upper and a lower bounds 
on the error respectively. For small deviations we can approximate the error using the second 
derivatives of M: 

A(pi V2 



(|>^ ^Ni{l - N,/N) 



(15) 

„ _ , (16) 

Rj ^Nj{l-Nj/N) 

To estimate the uncertainty induced by the specific bins choice, we preform all the analysis for 
a 1/2 unit redshift and LogiQ^L) binning and repeat for a 1/3 unit binning (where all bins widths 
are 1/3 unit, except the last two redshift bins which we cannot change because they contain too few 
data points). In the following, unless otherwise stated, we use the 1/2 unit binning for all further 
analysis. Clearly, the results with different binning are slightly different, but they are all within 
each other's error range. When we include the uncertainty induced by the binning the error ranges 
become only slightly wider. 



4 RESULTS 

So far, we have not assumed any functional form for the luminosity function or for the rate as 
our method does not require such an assumption. The "raw" results are depicted in Figures |2] and 
[3] Later on we will compare these step functions with models for the GRB rate that follow the 
SFR. However, in order to easily characterized the result we need to approximate our "raw" step 
functions with simple functional forms. Therefore, after obtaining the results in the form of step 
functions we approximate them in term of broken power laws: 

. L<L*, 
0(L)= (17) 
L>L* . 

Rgrb = Rgrb{0) ■{ ' ' (18) 

+ + Z>Zi, 
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Figure 2. The comoving rate and the observed events number redshift distributions 

Upper frame: The resulting step function and the best fit model: A broken power law with indices ni = 2.1, n2 = —1.4 and with a break at 
2; = 3.1 (solid line). The values for the models is 2.3 at 10 d.o.f. giving a rejection probability of 0.007. Lower frame: The number of detected 
bursts for each redshift bin and the rates expected from the model fitted above. The two upper boxes in each colunm represent the statistical en'or 
range. 

We obtain the parameters of the best fit functions by minimizing the values. 

The best fit functions are shown (together with the step functions) in Figures [2] and [3} The 
parameters of the best fit models are summarized in Table [1} To estimate the statistical errors 
involved in the parameter estimates we preformed a Monte-Carlo simulation. In this simulation 
we use the model with the best fit parameters to draw random sets of data (with same size as 
the real sample). We then carry out the analysis on this mock data sets and obtain new best fit 
parameters. Repeating this process many times we obtain a scatter of points in the parameters 
plane around the original best fit parameters. The central 68% and 95% ranges for each of the 
parameters separately, are also shown in Table [1] 

The luminosity function is well described by a broken power law, with a break at L* ~ 
10^^-^=^°-^[er5(/sec] and with indices a = 0.2lo;^ and (3 = IA^q^. The broken power law fit is 



very good, giving x = 0.63 for 4 d.o.f. This result agrees with previous studies (e.g. Daigne et al 
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Figure 3. The luminosity function and tlie observed and predicted luminosity distributions 

Upper frame: The resulting step function and the best fit model for the luminosity function: a broken power law with a break at L* = 10^^'^, a 
low luminosity index a = 0.2 and a high luminosity index /3 = 1.4. The value for this model is 0.63 at 4 d.o.f., giving a rejection probability 
of 0.04. Lower frame: The number of detected bursts for each luminosity bin and the rates expected from the model fitted above. The two upper 
boxes in each column, represents the statistical error range. 
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Table 1. Parameters results. The error ranges are 68% and 95% levels estimated using a Monte-Carlo simulations with 10000 sets, each of 101 data 
points. 
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2006, Guetta and Piran|[2007 ). The luminosity function cannot be fitted with a single power law, 



since such fit give high x^, rejecting such a model with high significance (98%). This contradict 



the results of Pelangeon et al. (2008 1 who studied the HETE-2 GRBs and found a consistency with 
a single power law luminosity function. The rate is described as well by a broken power law for 



1 + z, with a break at z = 3.1 



+0.6 
-0.8 



and indices of ni = 2.1^q g and n2 



The local event rate is po — ^■'i-o^[Gpc ^yr in agreement with previous studies e.g. 



Schmidt] ( fT999l )[po ~ 1.5] (see however |Schmidt| ( |200Tbl )[po ~ 0.15)]); [Guetta et al.| ( |2005| )[po 



0.5]; [Guetta and Delia Vallel ( |2007l )[po ~ 1.1]; [Liang etal] ( [2007] ) [po ~ 1.1]; [Pelangeon et al. 



(2008 )[po ^ 0.5]. The main factors determining the low redshift (current) event rate are the over- 
all GRB rate normalization, the low-redshift slope, the low end of the luminosity function slope 
and most important the position of the low luminosity cutoff. This low luminosity cutoff is essen- 
tial for any steep enough power law to prevent its divergence. This cutoff is critical to the question 
whether the model includes the low luminosity GRB population. The cutoff used in our estimates 
of pq'is L = lO^^erg/ sec. 



4.1 Consistency Checks 

Before we can accept the model,we turn now to check the validity of the assumption that the 
luminosity function and the rate are independent. . To do so, we preform a two dimensional 
Kolmogorov-Smirnov test (see [Fasano and Franceschini| ( |1987[ ), [Peacock[ ( [1983[ ), [Spergel et al. 



( |1987 )). Li this test the two dimensional data is compared to the modeled distribution for each of 
4 quadrant defined by axes crossing at each data point. The maximal difference is used to estimate 
the probability that the data is drawn from the distribution implied by the model. The probability 
that the models fits the data is displayed in the 2D K-S column in Table [4[ The test give high 
probability (96%), so we can accept the model and justify the underling assumption (of a redshift 
independent luminosity function). 

We carry out two other consistency checks. First, we compare the peak flux distribution ex- 
pected by our model with those observed by BATSE and by Swift. The results are shown in Fig- 
ures [4[ and |5j The KS-test resuhs, 83% for BATSE bursts and 17% for the fuU sample of Swift 
bursts (with or without redshift), indicate acceptable model. We note here, that the model is sig- 
nificantly rejected (KS-test < 10^) when compared to BATSE peak fluxes distribution using 
Piim = Q.2bph/ cm? / sec. However, when applying the effective detection threshold calculated by 



Band (2002) of Pum = 0.525p/i/cm^/sec, the model is accepted with high significance (83%). 
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Figure 5. Bursts count vs. peak flux for Swift bursts 

Cumulative bursts number distribution as a function of the peak flux p. Left: logarithmic scale showing A'^(< p) Right: log linear scale showing 
N{> p), which is used for the KS-test, giving a probability of 17%. 

Second, we compared the cumulative redshift distribution to the observed one and preformed a 
KS-test. Here, as well, the test gives a high probability for accepting the model (62%). 

To illustrate the distribution of bursts, we display the bursts in the L — z {Luminosity — 
Redshift) plane and in a rescaled plane in which the number of detected bursts (with or without 
redshift measurement) is proportional to the area (Fig. [6]). The rather uniform distribution on the 
rescaled plane is a visual demonstration of the validity of the assumption. 
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Figure 6. BursLs on the L-z plane 

The bursts distribution in the Luminosity-Redshift plane. Left: linear scale. Right: Axes are rescaled using R{z) and 4>(L), so that we have a 
uniform bursts-number density. The curved lines are contour lines for equal fluxes, for the values: p = 0.04, 0.4, 4, AOyph/ sec/ cm?], from top to 
bottom. The p = OA{ph/ sec/cm?] line is our detection threshold. The number density is uniform in the rescaled plane as expected and it drops 
to zero below the detection threshold. 



4.2 A Comparison With the Complete Swift Sample 

Our results are based only on the absorption and photometric determined redshifts of the Swift 
sample. Recently, Fynbo et al. ( 2009[ ) obtained emission lines redshift measurements for GRBs 
host galaxies that were not measured before. Most of those are at low redshifts. The growing 
number of emission lines redshifts in the range Q < z < 1 raises the question whether our model 
- based on a sample of redshifts measured from the afterglows - is consistent. We calculate, using 
our model, the expected redshift distribution of the entire bursts population and compare it to the 
number of observed bursts with a known redshift (measured using all methods). Clearly for any 
range of redshifts the number of known bursts with a given redshift range should not exceed the 
number predicted by the model. Fig |7] depicts this comparison for the cumulative number and 
for the counts number in redshift bins. Our model is consistent for all redshifts > 0.4. However 
there is an excess of low redshift bursts not predicted by our model. The discrepancy arise due 
to three bursts with z < 0.1 while less than one burst is predicted by the model. The redshift of 
the three bursts was determined using emission lines but no redshift was detected in this range 
using absorption lines. These three bursts yield a rate that is significantly higher than predicted 
by the model. This cannot be explained by a misidentification of the host galaxies because the 
probability for such an effect is too small ( Cobb and Bailyn|2008 1. These low redshift bursts have 
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Figure 7. Redshifts obtained by any method compared with the total number of bursts predicted by the model 

Upper frames: Left: The cumulative total fraction of bursts predicted by the model (dashed curve) and the detected cumulative fraction of bursts 
with redshift obtained by any method (solid steps) as a function of the redshift. Right: Zoom in on the low redshift part. Lower frame: The number 
density (per unit redshift) of predicted bursts (dashed curve) and an histogram of bursts with redshift obtained by any method (solid steps). Both 
are normalized by the total number of bursts. For all frames: The number of bursts with redshifts which had an at slew, i.e. the time elapsed form 
the trigger until the first optical observation is less than 300 sec (dashed dotted steps). 

low luminosities and they could not have been detected at a much higher redshift. We conclude 
that they possibly represent a low luminosity population that is different from the majority of the 



bursts (see a discussion in 6.2 1. 



4.3 The Redshift Distribution and Expectations for Future Missions 



A particularly interesting question is how many high redshift bursts are expected to be observed. 
Such bursts are of great interest as they may shed new light on the very early universe. Already now 
GRB090423 is amongst the most distant and hence the earliest objects observed so far. Figure [8] 
depicts the observed cumulative redshift distribution and the predicted one, for high-redshift bursts 
for Swift and for past and future missions: BATSE, EXIST ( |Bandet al.|2008| ) and SVOM ( |Schanne 
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Figure 8. The cumulative redshift distribution for Swift bursts and predictions for other missions 

The cumulative redshift distribution of Swift observed bursts and the model prediction (solid line). The KS-test gives a probability of 62%. Also 
presented are the predictions for some past and future missions: BATSE, SVOM and EXIST. Right frame is a zoom-in to the upper right part of the 
left frame. 



fraction of z > 


7 


8 


9 


10 


15 


20 


CGRO/BATSE 


0.4% 


0.2% 


0.1% 


0.05% 


0.01% 


0.002% 


Swift/BAT 


1.7% 


0.9% 


0.5% 


0.3% 


0.03% 


0.01% 


SVOM ECLAIRs/CXG 


2% 


1.2% 


0.7% 


0.4% 


0.06% 


0.02% 


EXIST/LET & HET 


4% 


3% 


2% 


1% 


0.2% 


0.06% 



Table 2. The high redshift fraction prediction for Swift and several other missions 



2008[|Gotz et al.||2009| ). The fraction of bursts with z > 7, 8, 9, 10, 15, 20 respectively is shown in 
Table H 

Detection of high redshift bursts is one of the main objectives of EXIST. The number of high 
redshift bursts expected to be detected by EXIST is presented in Table[3] During a five year mission 
EXIST will detect many high redshift {z > 10) (~ 30) bursts. In the best-fit model there is a good 
probability for EXIST to detect even a 2; > 20 burst during mission. This is, of course, provided 
that such early bursts exist and that the rate at z ^ 5 — 8 can be extrapolated to such high redshifts. 
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0.06 


+ 1.94 
-0.05 


0.01 


+0.79 

-0.007 



Table 3. High redshift detection rates prediction for EXIST and for SVOM. On average a redshift is obtained for only a third of the events. 



5 THE GRB RATE AND THE SFR 

The location of long GRBs in the star forming regions led to the expectation that GRB follow the 
SFR. We turn now to examine this hypothesis. Modeling the SFR throughout the measured redshift 
range (0 < z < 10) is a complicated task, involving observations in various wavelengths and 
various assumptions on observational proxies for the SFR as well as correcting due to obscuration, 



absorption and selection effects. In a classical work Madau (Madau et al. 1996, Madau 1998 



Madau et al. 1997 1 considered the SFR per comoving volume vs redshift and found a rise from 
present to 2; 1 — 2, and then a comparable decline to 2; ~ 5. This form of the SFR vs redshift has 
become known as the "Madau plot". A number of developments (e.g., dust corrections, submm 
results, new estimates of the SFR at low redshift) led to changes in the shape of the SFR. Following 
these developments Rowan-Robinson| ([T999 ) and Porciani and Madau ( 2001[ ) suggested that the 
SFR rises by a factor of 10 - 20 from 2; = to 2; ~ 1. The rate at higher redshifts {z > 2) was 
undecided at that time and models for the SFR at higher redshift included flat as well as rising and 
declining functions ( Porciani and Madaul|2001 ). Later, Hopkins and Beacom (2006) showed that 
at high redshift (z > 4) the SFR declines. Recently, several papers estimated the SFR using the 
new HST WFC3/IR camera dBouwens et al.|2009at |Oesch et al.|2009t [Bunker et al.|2009t [McLure 



et al. 


2009 


Yan et al. 


2009 



for ;z > 4 up to 2 ~ 8 — 9 . In the following we use Bouwens et al. (2009b) as representative of 
these new high-z SFRs. Despite all the observational advances, there are still different models of 
the SFR even at low (z < 1) and intermediate (1 < 2; < 3) redshift. On one hand, the widely used 



Hopkins and Beacom (2006) piecewise linear fit finds a factor 10 rise in the SFR from z = to 
2 = 1 and an almost constant rate from 2 = 1 to 2 = 4.5 (which follows by a steep decline at 
higher redshifts) .On the other hand [Bouwens et aL] ( |2009b[ ) use data from [Schiminovich] ( |2005[ ) 



and [Reddy and Steidel| ( |2009| ) for < 2; < 2, 2 < 2; < 3 respectively. Their data compilation 
can be fairly modeled (x^ = 10.8 for 9 d.o.f.) by a broken power law rising a factor of 20 from 
z = to z = 3.3, then declining as (1 + z)~^-^ for 3.3 < z. Note that a rise up to 2; = 2 — 3 is 



also suggested by Hopkins and Beacom] ( [2006] ) when fitting to a Cole (2001 ) form. This seemingly 
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mild disagreement between different SFR models have a crucial implication when comparing the 
GRB rate to the SFR. 



The same tests used in ^:4.1 can be used now to check the consistency of the data with a rate 
proportional to the SFR. An obvious problem is that SFR is not uniquely determined and hence 
we consider four different possible functions: Hopkins and Beacom ( 2006[ ) (denoted HE), SF2 



of IPorciani and Madau| ( |200T| ) (denoted PM SF2), |Rowan-Robinson| ( [T9991 ) (denoted R-R) and 



Bouwens et al. (2009b) (denoted B09). Note that HB and B09 decreases at large redshift, while 
PM SF2 and R-R stay constant. We compare, first, the SFRs to the binned rate we obtained by 
inverting the data (see Figure [9]) and we calculate the for both 1/3 and 1/2 binning. We find 
acceptable reduced values (see Table |4]) for all the functions. However, a comparison of the 
overall observed redshift and luminosity distributions with those predicted by models in which the 
GRB rate is fixed by a given SFR reveals that the 2D K-S or the K-S tests for the peak flux and 
the redshift distributions show inconsistency for the first three functions (HB, PM SF2, R-R). We 
find, however, consistency for the last one (B09). 

Next, we optimize the luminosity function for a given SFR. We take the GRB rate as known 



following a model of the SFR and obtain the best fit luminosity function by solving equation 12 
We now perform a 2D K-S test as well as K-S tests for the peak flux distribution and for the 
redshift distribution. The results of the statistical tests are shown in Table |4] Even though the fit 
improves still the first three SFR models fail the KS tests. The last SFR model (B09) is, of course, 
consistent. 

We attribute the consistency of the GRB rate with the B09 SFR model but not with the first 
three SFR models to the difference between the B09 SFR and the other three models in the range 
1 < z <?). While B09 keeps increasing in this range the first three are constant. This difference is 
crucial and understanding the SFR at this region is critical for the question whether the GRB rate 
follows the SFR or not. While there are differences at higher redshifts, the paucity of data leads to 
wide error range in that region allowing the GRB observations to be consistent with SFRs that are 
constant, decreasing or even increasing at large z. 

The comparison Table shows that the luminosity function parameters depends very weakly 
on the GRB rate model, a and po changes within their 95% error range, while L* and (3 stay 
almost unchanged. This illustrates the power of the method and the robustness of the luminosity 
function results, in particular for L* and (3. 
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redshift 

Figure 9. GRB event rate and several star formation rates 

Tile results for the rate, in 1/2 unit binning. Best fit for a broken power law - heavy black solid line. [Hopkins and BeacomH2006| SFR - red dashed 
line.'Bouwen s et al.|(2009b) SFR - cyan solid line. SF2 of |Porciani and Madau|{200T} - magenta dotted line. |Rowan-Robinson|(T999) SFR - blue 
dashed dotted line. 
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Table 4. Statistical tests for our models and for models following one of the SFRs considered (upper part) and for GRB rate from other studies 
(lower part). The consistent models are marked with a bold font. For each model we show the test with luminosity function from our results (first 
line) and with luminosity function that best fit observations after forcing the rate to follow the SFR (second line). 



5.1 Comparison with other works 



Our results are in agreement with most previous works on the GRB rate. Our results agree with 
Daigne et al.l (120061) who compared the Swift data with IPorciani and Madaul (1200 lb SFR models 
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and concluded that the rising model SF3 is the only one consistent (Our model is also consistent 
with SF3). We stress however, that the recent work on the SFR suggests that it decreases at high 
redshift and thus SF3 is, most likely, not a viable SFR model. Our results are also consistent with 



Guetta and Piran (2007 1 who find s that at 2; > 2.5 the GRB rate is significantly higher than the 



R-R or PM SF2 SFR. 

A comprehensive work studying the luminosity function and the rate of long GRBs was re- 
cently published ( Butler et al.|2009| ). While this work uses a different sample, adopts somewhat 
different assumptions and uses other methods, a comparison of the results can be very useful. We 
find a good agreement for the low luminosity power law index a = Q.22t^Q^^^{our results : 0.2]'^o;^) 
and for the break luminosity L* = 10^^'^^^° '^^(oMr results : io^2.5±o.2^ However we find differ- 
ent high luminosity power law index (3 = 2.9^^ (our results : 1.4to;6)- The differences may be 



explained by the fact that Butler et al. (2009) use a luminosity which is a time averaged whereas 
our luminosity is the peak luminosity. We find a nice agreement between the models for the GRB 
Butler etaL] ( |2009| ) find a rising rate for < 2; < 4 - first at slope 3.1 ± 0.7 for < 2; < 1 and 



rate. 



later at slope 1.4 ± 0.6 for 1 < 2; < 4. This is not very different from our results, recalling that the 
break at z = 1 was not a free parameter in their model. The decline slope for ^ > 4 is —2.9ll2^, 
but with the big uncertainties it also matches our model. It is thus not surprising that their model 
show consistency with the statistical tests (see Table|4]). Another useful comparison is with Kistler 



et al. (2009 1 who modeled the bias of the GRB rate with SFR in the range < z < A and used it 



together with the high-z bursts data to estimate the high-z SFR. The high-z GRB rate they found 
is roughly constant, in agreement with our results. 



6 SUMMARY AND CONCLUSIONS 

We find that the GRB rate increases up to redshift ~ 3 and it decreases at 2; > 3 (722 — —1-4), with 
68% confidence limits ranging from a steep decline (n2 — —2.4), to a positive incline (^2 — 1). 
The rate is compatible, of course, with a constant rate at higher redshifts. The model is consistent 
with all statistical tests and thus we can accept the basic assumption that the luminosity function 
dose not evolve with time. 

6.1 GRB Rate and the SFR 

The comparison between the SFR and the GRB rate seems to be inconclusive. This arises because 
of the differences between different estimates of the SFR at the intermediate redshift range 1 < 
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z <?). The rate we find is consistent with Bouwens et al. (2009b) (B09) (that follows Schiminovich 



2005) that describes a rising SFR from present up to 2; ^ 3. It is inconsistent with other SFRs 



( |Rowan-Robinson|[1999[ [Porciani and Madau||2001^ [Hopkins and Beacom||2006| ) that suggest a 
constant rate at this regime. 

At high redshifts the GRB data is sparse. The best fit result decreases slower than the most 
recent B09 SFR (or even the HB SFR) suggesting possibly higher GRB rate. However fast de- 
crease, like B09, cannot be ruled out while a flat or even slowly increasing rate are also consistent. 
A larger GRB rate at high redshift (compared to the SFR) can be explained within the framework 
of the massive stellar collapse model due to metallicity. Woosley and Heger ( |2006 ) suggests that 
GRB rate follow the low-metallicity part of the star formation. We cautiously note however that 



Fynbo et al. (2009) recently found that the optical afterglow spectroscopy sample is biased against 



measuring redshift at high metallicity environments, meaning that the result might be an artifact of 



a selection effect. Another clue can be taken from |Fruchter et aL] ( |2006| ) who found that the GRBs 
distribution within the host galaxies dose not follow the light distribution but rather some power 
(> 1) of the light distribution, i.e. higher GRB density in the denser star forming regions. This is in 
contrast to the core-collapse supernovae distribution that follows the light distribution, indicating 
the GRB progenitors are different from SN progenitors and hence their rate might be different. 



6.2 Low Luminosity Bursts 



We find an overall consistency when comparing our model and the full sample of GRBs with 
redshift. However, we also find three low redshift, low luminosity bursts (with emission lines red- 
shifts) that are not expected by the model prediction. The faintest burst in our sample GRB050724 
has a luminosity Li = lO^'^-^lerg /sec]. By applying our best fit model, we expect 0.9 bursts with 
luminosity L ^ Li, in the time span of our sample (4.5 years). Swifts weakest burst, GRB060218 
(Cusumano et al. 2006) with luminosity L2 = I0^'^'^[erg/ sec], is not in our sample because it 
has only emission lines redshift. Assuming that we can extrapolate the low end of our luminosity 
function we expect 2.5 ■ 10^® bursts with L ^ L2. Even when applying the 95% level of the pa- 
rameter a = 0.38, we expect no more than 5 ■ 10~^ bursts with L ^ L2. This implies that this burst 
represents a population of fainter GRBs, with much higher event rate, which cannot be directly ex- 
trapolated from the stronger GRB population (in agreement with e.g. jSoderberg et al.|2006[ [Cobb 



et al.[|2006t[Guetta and Delia Valle[|2007t [Liang et al.[[2007[ ). The emission lines redshifts sample 



includes two other low luminosity bursts: GRB051109 and GRB060505 with z = 0.08, 0.089 and 
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L = 10'^^-^[erg/sec], 10^^'^[erg/ sec]. The detection probability of such bursts, according to our 
model is < 2 ■ 10"'^, < 2 ■ 10~^ respectively. These three low luminosity bursts must belong to a 
different and a distinct group and we remove them from the analysis when checking for consis- 



tency in ^: 4.2 



6.3 High Redshift Bursts 

Higher redshift bursts are most interesting as they can provide clues on the very early universe. 
Extrapolating the rate for very high redshifts, we expect ~ 0.9 bursts with z > 8 and ~ 0.5 bursts 
with z > 9 (bursts with measured redshift), detected by Swift in the time span of our sample (4.5 
years). This is consistent with the observation of one z > 8 burst. Future missions like EXIST 



can find many more such bursts, even dozens of z > 10 bursts (see ^:4.3 1 provided that a simple 
extrapolation indeed hold to such high redshifts. 



6.4 The Local Event Rate 

The local event rate found is po = ^ ■•^-o'jlGpc'^yr'^], for bursts with L ^ 10^'^erg/sec . With 
one galaxy in lOOMpc^ this rate is equivalent to 1 event per galaxy per 10'' years! Taking into 
consideration the beaming factor of about 50 (see Guetta et al.|2005 ), the total events rate is about 
1 event per galaxy per 2 ■ 10^ years. This should be typical to our galaxy as a recent estimation 
to the SFR in our galaxy ( |Robitaille and Whitney1|2010| ) finds p = 0.68 - lA5MQ/yr, which is 
equivalent to the SFR in the local universe for a galaxy in a volume of lOOMpc^. One implication 
of this result will be on the possible association of GRBs with global extinctions of biological 
spices. These occurred on Earth a factor of 10 times less frequent, about once every 100 Myr. 
These two rates suggest that a typical Galactic GRB pointing to Earth does not cause a major 
extinction event. 
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APPENDIX A: THE LUMINOSITY 

In this work we used only data collected by Swift. The peak-flux of each burst, measured in Swiff 
BAT detectors band 15keV - 150keV ( Barthelmy et al.|2005 ). For most of the bursts we cannot es- 



timate the luminosity in the full 7-range (IkeV - lOMeV) as the spectral shape (the Band-function) 
( Band et aL]|1993 ) is poorly known. We still need however to have some quantitative measure of 



the luminosity that will be defined uniformly for all burst with a known redshift. To do so we use 
an average characteristic Band function, Ep^ak = 511A;eV^ (in source frame), a = —1, /3 = —2.25 
( Preece et al.|2000 ; Porciani and Madau|2001 ) to estimate luminosity, using the same parameters 



for all bursts, thus the estimated luminosity is proportional to the measured peak flux (p): 

= P^7cD{zYil + z)k{z)Cdet , (A.l) 

D{z) is the proper distance at redshift z, C^^^ is the fraction of the total 'j—ray luminosity detected 
in the detectors energy band for source at redshift = 0, 



!Z EN{E)dE 

Cdet = ^ 7^-^ , (A.2) 

k{z) is the k-correction for the given spectrum at redshift z, 

k{z) = .-^f'-;; — , (A.3) 

/a«,S:.::'v(sM£ 

where N{E) is the Band function, Emin = 15keV, Emm = 150keV. The values in this paper are 
the luminosity in the range [IkeV - lOMeV]. 
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Figure Al. Distributions of our fits parameters for tlie luminosity function (L*, a, (3) and for tlie rate (zi , ni , ?i2) wlien clioosing tlie 7-ray spectra 
(Band function) parameters randomly from their distributions. The shaded area contains the central 68% of the distribution, the dashed lines are the 
68% confidence range of our reported results. 



To study the dependence of our results on this approximation of an average Band function, we 
have preformed a simulation where the spectrum of the bursts is not universal but drawn randomly 
from the known distribution ( jPreece et al.||2000 ). Repeating the analysis in the paper 1000 times, 
each time calculating the luminosity of a burst using a band function with parameters randomly 



drawn - independently for each burst - from the distribution. Fig Al show the simulation results 
as an histogram for each of the parameters of the luminosity function and the rate fits. All the 
distributions concentrate within the 68% error range of the original results, demonstrating the 
robustness of the results and its insensitivity to the details of the distribution of spectral parameters. 



APPENDIX B: THE REDSHIFT DETECTION PROBABILITY 

We examine now the redshift detection probability as a function of the peak flux. We consider 
here two effects: first the probability to detect the GRB and second the probability to measure the 
redshift, for a given detected burst. 
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The probabilities to detect a GRB, or to measure a redshift, are a function of the burst energy, 
its duration and other factors, as described by (Band [2006 ). We restrict ourselves here only to 



the dependence of these probabilities in the peak flux, since this is the quantity we use for the 
analysis in this paper. 

Bl GRB Detection Probability 

The simplest model often used is of a sharp threshold: No detection with p < pum, but 100% 
detection of bursts with p ^ pum- For this model the value used at Swiff s main detection band 



[15 - 150]keV is pum = OAph/cm'^/sec. (see Guetta and Piran (2007), Gorosabel et al. (2004)) 



In the plot of the accumulated number of bursts as a function of log{p), a linear relation appears 
for fluxes that are low or medium but above the threshold mentioned above. Although we do not 
try to give a theoretical explanation for this result, we do not expect, however any strong deviation 
from that connection for lower fluxes, since our models predict a slow gradual smooth change in 
dN/dlog{p) and we can adopt this values at least as an order of magnitude estimators, to yield a 
continues threshold estimation. Assuming that the deviation from linearity for fluxes below pum is 
only due to a lowered detection sensitivity, we can extrapolate the predicted number of bursts for 
lower fluxes and extract the detection sensitivity by comparing the number of detected bursts to 



the predicted number. Figure Bl shows this: accumulated number of bursts vs. log peak flux, the 
linear fit and our fit. 



(l+c) , (1-c) 



erf{d ■ log{p/po) 0.2 ^ p 



p<0.2 

with the parameters found: 

c = 0.25, d = 10, loQioPo = -0.42 . 



Bl.l Redshift Measuring Probability 

Redshifts are measured only for a moderate fraction of the detected bursts. A few (5-10%) are too 
weak, some don't have a clear redshift signature and some are not measured because of lack of 
observational resources. 

When we consider the fraction of redshift-measured GRBs to the total number of GRBs de- 
tected, we see a trend of increase in this fraction with the measured peak number flux of photons 



in the detectors main band: 15keV — IbOkeV, p (see figure B2). We used a linear regression to 
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Figure Bl. The accumulated number of bursts as function of log(peak flux) 



approximate the relation: 

eMIO,{p) = a-log{p) + h, (B.2) 

Where , O^ijp) is the detection probability and 6z{p) is the probability that a redshift will be mea- 
sured. The parameters found in fit are: a = 0.080 ± 0.078 and b = 0.242 ± 0.054. with = 1-01 
at 5 degrees of freedom, giving a rejection probability of 0.038 . We expect the fraction of number 
of bursts with redshift Nz{p) to the overall number of bursts Nj(p), for any given p, to obey the 
relation: 

e,{p) N,{p) ■ ^ ■ ^ 

Figure |B2[ depicts the fraction of bursts with a measured redshift and the linear fit, for the 
absorption and photometry redshifts. Although the fit is acceptable the significant of the effect is 
just one standard deviation (a) - so with the current observations the option of no dependence of 
redshift detection probability with peak flux (a = 0) is still marginally consistent. 

Whereas the detection fraction can be fitted well with a linear model for the redshifts obtained 
using absorption lines, the case is different when considering redshifts obtained using emission 



lines. Figure B3 shows the measured redshift fraction and the linear fit for the emission lines 
redshifts. The emission lines redshifts are obtained preferably for high flux bursts, but very few 
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Figure B2. The fraction of bursts with a measured absorption lines and photometry redshift, as function of log(p) 
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Figure B3. The fraction of bursts with a measured emission Unes redshift, as function of log(p) 



are obtained for low and medium fluxes: 24% for logi^p > 1 and only 6.5% for logi^p ^ 1. This 
feature of the emission lines redshifts is associated with a strong bias toward lower redshifts as 
shown on Figure [T] These results supports our approach of selecting only the absorption lines and 
photometry redshifts as it make a sample which is less biased and easier to model. 
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52.54 


0.17 


1.32 


3.13 


1.85 


-1.38 
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52.53 


0.04 


1.32 
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52.53 


0.30 
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2.28 


-1.06 
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52.53 


0.17 


1.44 


3.11 


2.07 


-1.36 



Table Bl. Model's best-fit parameters when taking {■¥) and when not taking (-) each effect into account. 

The product of the above probabilities for burst detection and redshift measurement, gives the 
total probability that we detect a burst and measure its redshift 6^. Figure |B4| shows the function 
6z, which have the form: 

/ {alog{p) + + ^er/(d ■ log{p/p,)) 0.2 ^ p , 

I j9 < 0.2 , 

with 

a = 0.08, b = 0.242, c = 0.25, d = 10, logioPo = -0.42 . 



To check the effects of the burst detection probability and the redshift measurement probability, 
we repeated the process described in this paper, with all four options of taking or not taking 
into account each of the modified probabilities. The best fit results for the parameters are shown 



in Table Bl For both corrections, we find a non-negligible effect on the results, although the 



deviations induced on the parameters are in most cases within the statistical error ranges of the 
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analysis. When taking both effects into account, the models give somewhat better results in the 
various statistical tests. Therefor, the function 9z{p) that take both effects into account is used in 
this work. 



